Convergence conditions for iterative methods 
seeking multi-component solitary waves with 
prescribed quadratic conserved quantities 

T.I. Lakoba* 

Department of Mathematics and Statistics, 16 Colchester Ave., 
University of Vermont, Burhngton, VT 05401, USA 

September 25, 2009 



Abstract 

We obtain local (i.e., linearized) convergence conditions for iterative methods that seek 
solitary waves with prescribed values of quadratic conserved quantities of multi-component 
Hamiltonian nonlinear wave equations. These conditions extend the ones found for single- 
component solitary waves in [J. Yang and T.L Lakoba, Stud. Appl. Math. 120, 265-292 
(2008)]. We also show that, and why, these convergence conditions coincide with dynamical 
stability conditions for ground-state solitary waves. 



Keywords: Coupled nonlinear wave equations. Solitary waves. Iterative methods. 
PACS: 03.75.Lm, 05.45.Yv, 42.65.Tg, 47.35.Fg. 



lakobati@cems.uvm.edu, 1 (802) 656-2610 



1 



Solitary wave solutions of most nonlinear wave equations can be found only 
numerically. Recently, J. Yang and the present author obtained [Ij conditions 
under which an iterative numerical method can converge to stationary solitary 
waves of single-component Hamiltonian nonlinear wave equations. When this 
method, in what follows referred to as the imaginary-time evolution method 
(ITEM), converges, it provides one with a numerical approximation of a solitary 
waves with a prescribed value of a quadratic conserved quantity usually referred 
to either as power or the number of particles. However, many phenomena 
are described not by a single equation but by systems of coupled equations. 
Therefore, it is of interest to obtain conditions under which a mult i- component 
counterpart of the ITEM would be guaranteed to converge to, i.e., find, a multi- 
component solitary wave. We obtain such a condition in this work. Moreover, 
generalizing an observation made in [1], we show that the mult i- component 
ITEM converges only to those ground states of nonlinear wave equations which 
are dynamically stable, and explain why this is the case. 

1 Introduction and background 

For the single-equation case considered in Ref. [1], the power of a solitary wave is 

P= I u^(x)dx, (1.1) 



where u is the real-valued field of the solitary wave and x is the spatial coordinate. (Here 
and below, if the limits of the integration are not indicated, the integration is assumed to be 
over the entire spatial domain.) For example, if the time-dependent wave U{x, t) satisfies a 
Nonlinear Schrodinger-type (NLS-type) equation 

iUt + V^U + G{\U\^,x)U = C/(lxl ^ oo) ^ 0, (1.2) 

where V'^ is the Laplacian, then upon the substitution U{x,t) = e*^*u(x), where u is real, 
Eq. (|1.2p reduces to 

-fiu = L^^^u = 0, (1.3a) 

where 

L^^o)u = V\ + G{u\x)u. (1.3b) 

The parameter fi is referred to as the propagation constant of the solitary wave. A straight- 
forward calculation shows that the power P = J u^dx = J \ U\'^dx is conserved by the evo- 
lution equation (jl.2p . Thus, u can be parametrized either by P or by fi, so that one can 
write P = P{lj)- 

The method analyzed in [1] finds the solution u with a prescribed value of P by iterations: 

/^n = 1 ^ , (l-4a) 
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Un+1 -Un = [l^^^^Un " /^n^^n) Ar , (1.4b) 

/ P 

Un+l = Un+l\ jz r. (1.4c) 

where Ar > is an auxiliary parameter, and a positive definite operator A'^ can be conve- 
niently chosen in the form mimicking the linear constant-coefficient part of L^^^ in p.3ap : 

iV = c-V2, c>0. (1.5) 

The purposes of Ar and N^^ in ()1.4bp will be clarified shortly. The inner product is defined 



as 



(/(x),5(x)) ^ j{f{^)Yg{^)d^. (1.6) 

Methods similar to (II. 4p had also been considered for finding solitary waves in the past 
(see, e.g., [2]-|5]). However, it was in ^ where the convergence conditions of the specific 
version, (jl.4|) . of the ITEM were found in terms of the properties of the linearized operator 
of Eq. (|1.3a|) . Namely, Eqs. (|1.4p are linearized by a substitution 

Un = U + Un, \Un\<.U, (1.7) 

which results in [T] 

AT-'^r- A r~ - T~ {N-^u,Lun) 

Un+l -Un = N LUn^T, LUn = LUn 7T— j ^ U . (1.8) 

{N~'-u, u) 

Here L is the linearized operator of L^^^ in Eq. (|1.3ap . For example, for L^^'^^ given by 
(fObj) . 

L = - ^ + G('u^ x) + 2u^G^2 (u^ x) . (1.9) 

The ITEM (jl.4p converges if the eigenvalues of the right-hand side of its linearization (jl.Sp 
are located between —2 and 0. The first of these conditions implies that 

Ar < (1.10) 

where Amin is the most negative eigenvalue of N~^C (In practice, the maximum Armax can 
be easily found by just a few trials, and then the value Ar leading to an optimal convergence 
rate is usually somewhat smaller than Armax! see Fig. 3 and Eq. (47) in [Ij.) The second 
condition can be shown [1| to be related to the properties of the original equation (jl.3ap 
and its linearized operator L as follows First, let us denote 

p{L) = the number of positive eigenvalues of L counting their multiplicity (1-11) 

(and similarly for any other operator or matrix). Next, assume that we are considering the 
generic case whereby the null space of L does not contain any functions other than those of 
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the modes du/dx^^^ which are associated with translational invariance of the sohtary wave 
along x^''). Then, under condition (fLTOll . the ITEM (O]) converges provided that either 



p{L) = and ^ 



(1.12a) 



or 



p{L) = 1 and > . 



(1.12b) 



In all the other cases algorithm (|1.4|) diverges [9]. Remarkably, as pointed out in [T], these 
convergence conditions are the same as the stability conditions of nodeless solitary waves in 
the NLS-type evolution equation (jl.2p [8j. In other words, the ITEM ()1.4p converges only 
to those nodeless solitary waves of (|1.2p that are dynamically linearly stable. (Let us note, 
in passing, that iterative methods that are guaranteed to converge to any solitary wave, 
stable or unstable, were first proposed in [3] and later developed in [TO].) 

The purpose of using operator A^~^ in ()1.4bp is to considerably reduce the magnitude 
of the most negative eigenvalue of operator N^^C compared to that of L [3l[Tl[Tn]. This is 
analogous to preconditioning a poorly conditioned matrix A when solving a linear equation 
Ay = b by an iterative method; this can considerably improve the convergence rate of 
the iterations (see, e.g., [7], Lecture 40). In regards to implementing N, note that it is 
a differential operator with constant coefficients and hence has a simple representation in 
the Fourier space. Therefore, N'^L^^^'^Un and Un are easily computed using the direct 
and inverse Fast Fourier Transforms, which are available as built-in commands in all major 
computing software. 

Note that the propagation constant fi in the ITEM (|1.4p is not prescribed but computed 
iteratively. Numerical methods for finding solutions of Eq. ()1.3p with a specified value of fi 
rather than with a specified value of power (jl.ip have also been considered in quite a few 
studies (see, e.g., references in |llj). and we will not consider them in this work. 

In this paper, we derive convergence conditions of a generalization of the ITEM (11. 4p 
for multi-component solitary waves. Such a generalization was proposed in [TTj (but see 
also Section 4.3 in [TO]), and its algorithm is presented in Section 2. Note that in [11], we 
also proposed another method that finds solitary waves with the same conserved quantities 
as the ITEM but converges faster; moreover, it converges much faster than the ITEM when 
the latter converges slowly. This method is a modified form of the well-known Conjugate 
Gradient method (CGM; see, e.g., [7], Lecture 38), and its algorithm is given in Appendix 
for the reader's convenience. In [TTj we showed that this modified CGM has the same 
convergence conditions as the ITEM. Therefore, these convergence conditions, which gener- 
alize conditions (jl.l2p . apply to both the ITEM and modified CGM. The derivation of these 
conditions is the main result of this paper and is presented in Section 3. This derivation 
is based on the idea of Ref. [12], where it was used to establish stability conditions for a 
certain class of multi-component solitary waves. In fact, our convergence conditions of the 
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ITEM and modified CGM turn out to be the same as the stabihty conditions derived in [12j. 
This relation between the two sets of conditions generalizes a similar fact pointed out in [l] 
for single-component equations, and in Section 5 we explain under what circumstances such 
a coincidence of the convergence and stability conditions occurs. In Section 4 we provide a 
geometrical argument that facilitates intuitive interpretation of the convergence conditions 
derived in Section 3 for the special case where the solitary wave has two quadratic con- 
served quantities. In Section 6 we summarize the results of this work. Let us emphasize 
that numerical examples involving the multi-component ITEM and CGM are not a subject 
of this analytical study; the interested reader can find such examples in Ref. [TT] . 



2 ITEM algorithm for mult i- component solitary waves 

Let us begin with an example that will motivate introduction of some new notations. The 
following system describes pulse evolution in a two-core nonlinear directional fiber coupler 
where each core supports two orthogonal polarizations of light |13j : 

+ ^(1) + (^|^(l)|2 + ^|f;(2)|2^ f^d) + ^(3) ^ 
) + + ( |[/(2)|2 + ^|C/(1)|2) [/(2) + f;(4) ^ q 



+ 



2 + ^|c/(4)|2j um^uw^Q 

2 + ^|C/(3)|2\ [/(4)+^(2) ^0 



(2.1) 



Here {U^^\ U^'^^) and {U^^\ U^^^) are the pairs of orthogonal polarizations in the two cores. 
The quadratic quantities conserved by these equations and generalizing (jl.ip are: 

' p(l) + p(3) \ / 1 1 \ 

p(2) + p(4) 1 - I 1 1 j 



Q 



p(2) 
p(3) 
p(4) 



(2.2) 



where P^'^) = {{U^^^)*, U^^^) . Upon the substitution 



\ 



(2.3) 



/ u^^\x) 
u(^\x) 
u^''\x) 
mW(x) j 

where u^^^ can (for the purpose of this example) be chosen to be real-valued, system ()2.1 
reduces to: 



V 



^il'i + ((n(2))2 + K(tx«)2)n(2) + n 

.,(4) , /'^,(4)^2 , ^^,(3)^2^ „(4) , ,,(2) 



(4) 
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(2.4) 
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where jl = (/^^^\ //^^^)"^. Then the powers of the individual components are 

= (n('=),-uW). (2.5) 

We now generahze this example to an 5-component system possessing an s-component 
vector of conserved quantities Q, so that the kth component of Q is: 

s 

Q{k) =Y^q(kl)p{l)^ fe = l,...,s<5, / = !,..., 5, (2.6) 

1=1 

where the solitary wave is u = [u^^\ . . . jU^"^))"^. As illustrated in the above example with 
S = 4 and s = 2, the number of conserved quantities can be less than the number of 
the components of the solitary wave: s < S. Other examples where the situation s < S 
takes place include: a system of NLS-type equations coupled coherently via phase-sensitive 
nonlinear terms (as opposed to linear ones as in (j2.ip ): the well-known system of three 
waves interacting via quadratic nonlinearity [14:\ I15j : and any system of coupled carrier- 
wave (also known as long-wave, or Korteweg-de Vries-type (KdV-type)) equations, as we 
will explain at the end of Section 5. To emphasize the possibility of having s < 5, we use 
a different vector notation for Q than for u. Now, the matrix (g^^'^) in (|2.6p is assumed 
to be in reduced echelon form (see any textbook on undergraduate Linear Algebra) and, in 
addition, its columns are arranged so that 

In ([2:2]) . matrix [q^^^^] is the 2 X 4 matrix. 

The multi-component generalization of Eqs. (|1.3a|) and p.4ap is 

L(o)u = l(°°)u - U {N-^U, U)-^ (N-i^, L(°°V) = , (2.8a) 

U^^-^. (2.8b) 
Ou 

For example, in (lOll . L(°°) u is the first term (the 4x1 vector), and U is the first factor of 
the second term (the 4x2 matrix) on the left-hand side. The SxS matrix N is a self-adjoint 
positive definite operator. For optimal preconditioning, its differential part should mimic 
the highest derivative in the linear part of L'-^^^. For example, for the L*-^*^^ in ()2.4p . N can 
be chosen as a diagonal matrix with its diagonal entries of the form (II. 5p . 
The multi-component version of the ITEM (jl.4p is: 

u„+i - = (l(°°)u„ - Un {n-^Un, Un)-^ {^-'Un, L^^^V^)) Ar , (2.9a) 



(fc) _ .(fc) 



. fc = l,...,.<5, (2.9b) 

\ ^n+l 
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where 

^n+l — \^7i+l; ""n+l/ ' K — i, . . . , .5 . 

Note that, by (j2.7|) . the numerator of the fraction under the square root equals P^^\ Let 
us emphasize that if s < 5", then ()2.9bp specifies only that the s components of Q have 
their prescribed values but does not impose any other conditions on the powers, PW, of 
the individual components of the solitary wave. 

As we noted in Section 1, in [llj we proposed a modified CGM that converges under 
the same conditions that we will establish for the ITEM (|2.9p . but faster. Moreover, it 
converges much faster when the ITEM converges slowly. Its algorithm, however, is somewhat 
less transparent than (j2.9p . and therefore we state it in Appendix. Let us reiterate: The 
convergence analysis that we will present in Section 3 applies both to the ITEM and the 
modified CGM. We advocate using the latter method when the slow convergence of the 
ITEM justifies spending a little extra effort on programming the algorithm of the CGM. 

As in the single-component case, we perform convergence analysis of the ITEM p.9p 
using linearization analogous to (II. 7p . A tedious but straightforward calculation shows that 
the linearized operator of the right-hand side of Eq. (j2.9ap is: 

n-^Cvin = (La„ - u {n-^u, uy^ (n^^z^, Lii„) ) . (2.10) 

Here L is the linearized operator of L^*^) in (j2.8ap obtained when the last term in that 
equation is replaced by lA fl; compare with (j2.4p . (Although fl is not prescribed but instead is 
iteratively computed within the method, its exact value can still be used in the convergence 
analysis.) For Hamiltonian wave equations, L is self-adjoint. Next, the conservation of Q 
implies the orthogonality relation 

{U,u^) = 6. (2.11) 

Taking the inner product between U and (j2.9ap . one can show that Eq. ()2.9bp does not 
change the linearization of p.9ap : the role of (j2.9bp is to guarantee that the s components of 
vector Q equal their prescribed values exactly rather than in the linear approximation. Thus, 
the operator in ()2.10p is the linearized operator of the multi-component ITEM. Operator C 
is easily verified [J to be self-adjoint on the space of functions satisfying the orthogonality 
relation (12. lip . However, N^^£ is not self-adjoint. To cast the linearized ITEM ()2.9p into a 
form involving only self-adjoint operators, which is more convenient to analyze than (j2.10p . 
we use the following change of variables: 

v„ = n1/2u„, V = N"^/2z^, /C = N~^/2£n-i/2^ K = N'^/^lN-^/^ (2.12) 

Then the linearized ITEM (j2.10p and the orthogonality relation (|2.1ip become: 

v„+i - v„ = /Cv„Ar, /Cv„ = Kv„ - V (V, V)"^ (V, Kv„) , (2.13) 

(V,vj = 0. (2.14) 
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In what follows we will analyze the transformed form (j2.13p of the linearized ITEM, 
because it involves operator /C that is self-adjoint on the space of functions satisfying the 
orthogonality relation (j2.14p . Therefore, the evolution of the iteration error v„ is completely 
determined by the eigenvalues of /C. For convergence of the ITEM, these eigenvalues must 
lie between — 2/Ar and 0. The first of these conditions is achieved by adjusting Ar, whereas 
the second condition is analyzed in the next Section, where we will establish its connection 
to the number of positive eigenvalues of L. It should be pointed out that by Sylvester's law 
of inertia (see, e.g., [H]), the numbers of positive and zero eigenvalues of K and L are the 
same. Therefore, we will refer everywhere to those eigenvalues of L, whereas in the analysis 
of ()2.13|) it is the eigenvalues of K that are involved directly. 

Finally, a note is in order about the effect of zero eigenvalues of L. As in [1], we assume 
the generic situation whereby the null space of L does not contain any functions other 
than those of the modes du/dx^^^ which are associated with translational invariance of the 
solitary wave along coordinate x^'^^ As was shown in [Ij and [llj for the single-component 
ITEM and CGM, such modes lead only to a slight shift of the solitary wave along the 
respective coordinates, but do not otherwise affect convergence of the iterative method. 
The same proofs carry over directly to the multi-component case. Thus, in what follows we 
will focus on nonzero eigenvalues of L. 

3 Stability criterion for mult i- component iterative methods 

This Section contains the main result of this work, Eqs. (|3.ip . which are derived using 
a combination of analyses of Refs. [I] and [12]; see [TTj. Namely, we will show that the 
operator JC is negative definite on the space of functions satisfying ()2.14p provided that the 
Jacobian matrix 

--— = ^, IS nonsmgular (3.1a) 

and that 

p(L)=p(^), (3.1b) 




where the notation p is defined in (|l.lip . These conditions generalize conditions (jl.l2p 
for the multi-component ITEM (j2.9p and modified CGM ()A.ip . As explained at the end of 
Section 2, under these conditions the ITEM can be guaranteed to converge by choosing Ar 
to be sufficiently small (in practice, Ar = 0(1) [llllO]). The modified CGM is guaranteed 
to converge provided that (|3.ip hold. 

Let ^ be an eigenfunction of /C and $ be an eigenfunction of K: 

/C^- = A^', K$ = A$ . (3.2) 

Taking the inner product of V with the first equation in (j3.2p and using the definition of 
/C, one sees that eigenfunctions ^ with A 7^ satisfy the orthogonality relation (|2.14p . 
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However, the eigenfunction ^ with A = does not, in general, satisfy that relation, as we 
will see later on. Let us now expand ^ and V over the set of <I>'s: 



= ^am'^m{^) + / a(A)$(A,x)dA, 

J cont 



(3.3) 



Here the two terms in each expansion correspond to the contributions of the discrete and 
continuous spectra of K, and a's are scalars and B's are s x 1 vectors: 

Bm = {V,^m), B{X) = {V,^X)). (3.4) 

Here and below we do not indicate the dependence of ^ and ^' on x since it is always 
implied. Let us also denote an s x 1 vector 

^ = (V, V)-i (V, KM/(A)) . (3.5) 

From the ([33]) one finds: 

am = BlH/{Xm-A), a{X) = B^{X)H/{X-A). (3.6) 

Substitution of these equations into the orthogonality relation ()2.14p . which is to be satisfied 
by ^'(A) for A 7^ 0, yields: 

Before continuing, we need to point out one important feature of the eigenvalue problem 
for ^'(A), which can be restated as 

- A"^ = VH . (3.8) 

Namely, the vector H in it is arbitrary, and, therefore, by specifying different H^s one 
obtains different ^''s for a given A. To verify this, one only needs to substitute K^' from 
(j3.8p into (j3.5p . This observation about H being arbitrary allows one to reformulate the 
problem R{A)H = as follows: 

Analyze under what conditions matrix R{A) can be singular for A > 0. (3.9) 

If we find that R{A) can be singular only for A < 0, this would imply that /C is negative 
definite (modulo the remark made at the end of Section 2) and hence the multi-component 
ITEM and modified CGM would converge. Again, note that in arriving at formulation 
(j3.9p . we have relied on the arbitrariness of H. Indeed, if H had not been arbitrary but 
instead determined by ^'(A), then the first equation in (|3.7p would not have been equivalent 
to (j3.9p . since even though -R(A) could have been singular, the particular H might not have 
necesarily been its eigenvector. 
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To address question ()3.9p . we study the eigenvalue problem for the s x s real and 
symmetric matrix R: 

R(p = -f{A)0. (3.10) 

Its eigenvalues can be found from the Rayleigh quotient: 

7(A) = ^^i?(^/<^^(^. (3.11) 

From ()3.1ip . ()3.7p . ()3.4p and the completeness of the set of <I>'s one has: 

7(A ^ +00) = - {V^, Vip) I {'fvK) -0 . (3.12) 

Here we have used the fact that rank of V is s, since V is constructed from s independent 
components of Q] hence V(p 7^ 0. Similarly, one verifies that all the eigenvalues of R satisfy 

d-f/dk < when A / A^. (3.13) 

As A — > Xm, the matrix (Am — A)i? — > BmB^. The latter is a rank-one matrix, and 
hence only one of its eigenvalues is nonzero. Therefore, at A = Am, at most (see below) 
one eigenvalue of R has a simple pole singularity and the other eigenvalues are continuous 
functions of A. 

The facts stated in the previous paragraph allow one to specify when R can be singular 
(i.e., one of 7(A) = 0) for A > 0. We will do so for generic cases first and at the end will 
consider the missed special cases. One should consider three possibilities: 

(i) s = p{K), (ii) s>p{K), (iii) s<p{K). (3.14) 

Qualitatively different situations for possibilities (i) and (ii) are examplified by Fig. [T]^a-d). 
From panels (a,c) in this Figure one can see that 7 = does not occur for A > provided 
that p(K) = piR{0)). From Fig. TO^,d) it also follows that if p(K) > p{R{0)), then there 
is always a 7 = for some A > 0. By inspection, one can convince oneself that these 
statements are true is the general case (i.e., for any s and p(K)) for possibilities (i) and 
(ii). One can also see that the situation where p(K) < p{R(0)) cannot occur. Indeed, by 
(j3.12p . all 7(A +00) < 0, and they can become positive only at Am's. Hence the number 
of positive eigenvalues of -R(O) cannot exceed the number of positive Am's, which is p(K). 
Finally, in regards to possibility (iii), one can easily see (Fig. [T]^e)) that there should always 
be a 7 = for some A > 0. To summarize, R does not become singular for positive A only 
when p{K) =p{R{0)). 

Thus, to arrive at conditions (j3.ip . we need to relate R{0) with dQ/dfl. First, in analogy 
to the first equation in ()3.7p . 

i?(0)i? = (V, ^(0)) . (3.15) 

As we noted after Eq. (j3.2p . the right-hand side of (j3.15p does not, in general, vanish. We 
will now find ^'(0). From the definition of /C, 

K^(0) = V#. (3.16) 
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Figure 1: Various cases showing the eigenvalues of matrix i?(A), as explained in the text. 
The shaded oval along the vertical axis indicates that the curves which are depicted as 
crossing that axis at positive values can actually cross it anywhere (i.e., also at negative 
values), (a) s = p{K) = 2; (b) s = p{K) = 2; (c) s = 3, p(K) = 2; (d) s = 3, p{K) = 2; 
(e) s = 2, p(K) = 3; (f ) s = p(K) = 3; (g) s = p(K) = 3; (h) s = 2, p(K) = 3. 



11 



On the other hand, consider Eq. (j2.8ap written for the transformed operator K^^'^-'v = 
N~^/^L(°''^u and note that the last term in that equation is Vjl (see also (j2.4p ). Differen- 
tiation of this equation with respect to fj.^''^ yields: 

K^ = Ve('=), k = l,...,s, (3.17a) 

where v = N-'^/^u and e^'^) is the s x 1 vector whose A;th entry is 1 and the other entries are 
zero. Combining Eqs. (|3.17ap for all k yields 

k| = V, (3.17b) 

where dv/dfl is an 5 x s matrix. Now, comparison of (13TB and (I3.17bll shows that 

^r(O) = dv/dfl + 5^'^ dv/dx^^^ , (3.18) 

where g^^^ are arbitrary constants, and we have used our assumption that the null space 
of L (and hence of K) can consist only of translational-invariance eigenmodes. All such 
eigenmodes are orthogonal to V, as can be seen by considering their inner products with 
()3.17bp . Then the substitution of (j3.18p into (|3.15p and recalling that V = 6Q/6v (see 
(IZiSbl) and ([21^ ) yields R{0)H = dQ/dfl H. Given the arbitrariness of H (see the text 
after ()3.8p ). this implies that 

R{0) = dQ/dfl. (3.19) 

This fact along with the summary sentence found before Eq. (j3.15p entails condition ()3.1bp 
in the generic case. 

Let us now consider special cases that we have glossed over. First, suppose one of the 
terms in the discrete sum in (j3.7p with a Am > is a zero matrix. This can occur only 
when B = 0. Then by the first equation in (|3.4p , the corresponding eigenfunction of K 
satisfies the orthogonality condition p.l4p and, by ()3.8p and ()3.5p . is also an eigenfunction 
of /C with the eigenvalue A = > 0. Thus, even though all the eigenvalues 7 of i?(A) 
are continuous at A = A^ and do not change their signs (see Fig. [U^f), where Am = A2), 
operator /C still has a positive eigenvalue A = Am- Note that since fewer than p(K) of the 
eigenvalues 7 change sign as A decreases from +00 to 0, then p{R(0)) < p(K) and hence by 
(j3.19p . this special cases falls under the generic condition ()3.1bp . 

Second, suppose that two positive eigenvalues of the self-adjoint operator K are the 
same. Since the corresponding eigenfunctions are linearly independent, so are the eigenvec- 
tors of R whose eigenvalues will have a pole singularity at A = Am- Figures [T]J^g,h) illustrate 
two qualitatively different situations that can occur in this case. As one can see, condition 
(|3.1bp still determines whether any of the 7(A) 's can vanish for A > 0. 

Third, suppose dQ/dfl = R{0) is singular. This means that there is an eigenfunction 
^'(0) of /C that satisfies the orthogonality condition (j2.14p . As we point our below, this 



12 



may prevent the ITEM and modified CGM from converging. Hence we impose condition 
(j3.1ap . Thus, we have established both conditions ()3.ip as being necessary and sufficient 
to guarantee that the iterative methods (j2.9p (with (jl.lOp being satisfied) and (jA.ip will 
converge for any initial guess uq that is sufficiently close to the solitary wave being sought. 

Let us now discuss how these methods may behave if either of these conditions is vio- 
lated. If ()3.1bp does not hold, then the ITEM is guaranteed to diverge for a generic initial 
condition, since the iteration error u„ will contain a component of the eigenmode that will 
increase by a factor (1 + AAr) > 1 at each iteration. On the other hand, if (j3.1ap is vio- 
lated, then the iteration error will contain an eigenmode with A = 0, which will remain the 
same at each iteration. Hence the ITEM will settle near (u + small constant • ^'(0)), and 
the norm of the iteration error will not be able to reach an arbitrarily low prescribed error. 
Thus, the linearized convergence analysis predicts that the method will "stall" at a higher 
error, but will not diverge. (Taking into account terms nonlinear in Un (see (jl.7p ) may yield 
the information of whether the method actually converges or diverges. However, such a 
nonlinear analysis is of limited practical use since, even if the ITEM is found to eventually 
converge, it would do so very slowly in this case.) 

As for the modified CGM, it is not bound to diverge if IC is not negative definite. 
However, it may do so when either of conditions (j3.ip is violated. The mechanism of this 
divergence would be the vanishing of the denominator in Eqs. (lA.lb .e) of the algorithm. 
Such a divergence may perhaps be avoided by choosing a different (but still generic) initial 
guess uq. 

Finally, let us present an example where one can predict convergence of the ITEM and 
modified CGM without computing the eigenvalues of L and dQ/dfl. This example is a 
straightforward extension of Corollary 1 in [Ij to the multi-component case. Consider a 
system of incoherently coupled NLS-type equations, generalizing (|1.2p : 

^[/W + V^C/W + GC^)^!;/^!^,... ,|t/(^)p,x) C/W = 0. (3.20) 

Its solitary wave is sought in the form U^''\x,t) = u^^\x.) exp{ifi^^H) with u^^^ being real. 
Note that in this case, s = S,Q= (P(i) , . . . , P(^)) , and also 

L = + G, (g)('=') = 2n('=)n« dG^''^ /d{u^'^f . (3.21a) 

In this case the operator L^''^ is diagonal with entries 

(lW^W^.^W + V^ + G^'^). (3.21b) 

Suppose that at least one of the components of u has at least one node and Q is positive 
(semi)definite. Then p(L) > S, and hence the iterative methods (j2.9p and (jA.ip diverge. 

On the other hand, suppose that all components of u are nodeless and G is negative 
(semi)definite. Then p(L) = 0, and hence the iterative methods converge. The proof of 
both statements repeats that of Corollary 1 in [IJ and hence is not given here. 
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4 Geometric interpretation of p{dQ/djl) when s = 2 

In the case of a single equation, the ITEM's convergence condition dP/d^i > in (|1.12bp 
has a simple geometric interpretation: the curve -P(^) must have a positive slope. (In other 
words, solitary waves for which p{L) = 1 and curve -P(/i) has a negative slope cannot be 
obtained by the ITEM; however, they may be obtained by other iterative methods [I].) The 
convergence condition in the multi-component case, Eq. ()3.1bp . is not as straightforward to 
visualize. In this section we will give a geometric interpretation of this condition for the case 
s = 2, i.e. when the solitary wave has two quadratic conserved quantities Q*-^^ and Q^'^\ 
(Note that the solitary wave in this case can have more than two components: examples 
include the three- wave system |14tll5j and Eqs. ()2.ip .) While in the s = 1 case the Jacobian 
dQ/dfl = dP/dfi can be either positive or nonpositive (i.e., there are two possibilities), in 
the s = 2 there are three possibilities: when dQ/dfl has two, one, or no positive eigenvalues. 
Thus, below we will give a geometric interpretation of these three situations. Interestingly, 
this interpretation makes reference of a single curve (see Eqs. (|4.5p below) — an intersection 
line of surfaces Q «(^«,/x(2)) and (^(1)^^(2)) shifted in a certain manner. 
For brevity, let us denote 





an 


Ol2 \ 






022 / 



(4.1) 



Note that dQ/dfl = R{0) is a symmetric matrix (see p.7p ). and so 021 = ai2- From the 
quadratic equation satisfied by its eigenvalues one can see that: 

det{dQ/dfi) > and 011,022 > p{dQ/dfl) = 2; (4.2a) 

det(5Q/5/2) < p{dQ/dfl) = 1; (4.2b) 

det(5Q/a/2) > and 011,022 < ^ p{dQ/dfl) = 0. (4.2c) 

We have used the symmetry of dQ/dfl to infer that a definite sign of (on +022) in (j4.2b .cl 
implies the corresponding sign for both these diagonal entries individually. 

Let us consider two surfaces Q^^^ fi^'^^') and Q*-^^ (^'■^^ /U*^^^) and the normal vectors 
to them at a given point fi^"^^) : 

= [ofci, Ofc2, -1] , k = l,2, (4.3) 

where these vectors are chosen to point downward. If these surfaces are vertically shifted 
so as to have the same height at point (/u^^\ /i*-^-*) , then the cross-product of the normal 
vectors defines the direction of the intersection line of such shifted surfaces at this point. 
The vertical component of this cross-product is det{dQ/djl): 



rii X rir 



ail ai2 -1 
021 022 -1 



(4.4) 
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Thus, according to the above definition, the intersection hne, of the shifted surfaces 
Q^^^ fi^^^) and Q'-^-' (^^"^\ /i*-^-*) points in the same vertical direction (i.e., up or down) 
as X ng. Let us also note that (a^i, 0^2) are the projections of n^j. onto the axes 11^^^ and 
/i(^). With these observations, conditions ()4.2p are restated as: 



i points up and 

projection of rij^. on respective axis fi^''^ is positive 

i points down 
I points up and 

projection of on respective axis fi^'^^ is negative 



pidQ/dp) = 2; (4.5a) 



p{dQ/dfl) = 
p{dQ/dn) 



0. 



(4.5b) 
(4.5c) 



Finally, let us note that conditions (j4.5p can be restated solely in terms of the two- 
component vectors = [aki^o,k2\ obtained by projection of on the horizontal plane: 



angle between ai and a2 is < 180° and 
projection of on respective axis /U*^^-* is positive 



» p{dQ/dfl) = 2; (4.6a) 

> p{dQ/dfl) = 1; (4.6b) 

^ p{dQ/dfI) = 0, (4.6c) 
where the angle is measured from oi to 02 in the counterclockwise direction. 



angle between oi and a2 is > 180° 

angle between ai and 02 is < 180° and 
projection of Ok on respective axis fi^^'^ is negative 



5 Connection between convergence and dynamical stability 

First, we observe that the conditions (j3.ip under which the ITEM (j2.9p and the modified 
CGM ()A.ip are guaranteed to converge [18] are the same under which the solitary wave of 
the incoherently coupled NLS-type equations p.20p with all nodeless components is linearly 
stable [12]. (More precisely, the system analyzed in [12J had G^^^ = Y^i=iCrki\U^'-^'^ , but 
the results of that paper are straightforwardly extended to apply to (j3.20p .) Thus, the 
nodeless solutions of (|3.20p found by the iterative methods of this paper are guaranteed 
to be dynamically linearly stable. This is an extension of the result found in [IJ for a 
single-component Eq. (jl.2p . 

Let us now explain why this close connection between the convergence and stability 
takes place. Our explanation applies both to single- and multi-component equations. In 
regards to the convergence conditions, recall that the iteration methods converge when 
the operator C in (j2.10p is negative definite on the space of functions V satisfying the 
orthogonality relation ()2.1ip . Note that on this space, (ipjCip) = (ipjljip), and therefore, 
the negative definitenesses of C and L are equivalent under (12. lip . 

Let us now turn to the stability conditions. The details are slightly different for envelope 
and carrier solitary waves, so we begin with the former using Eqs. (j2.ip as an example 
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whenever needed. Seeking the shghtly perturbed sohtary wave in the form similar to ()2.3p 
where now all u^^^ are replaced by 

+ + inf))e^*, ng),(x) € M, (5.1) 

one obtains (see, e.g., [E]): 

Lur = Xuj, Ljuj = -Xun . (5.2) 

In the case of Eqs. ([23]) or UjTlOli . Lj = L(°) (see dZHa])), but in general (e.g., for the three- 
wave system |14[I15]) this is not so. Also, in the case where s < 5, it is convenient, although 
not critical, to write the /2-term in L/U/ as, e.g., for (|2.4p : diag(^fi^^\ fi^'^^) u/; this 

makes Lj explicitly self-adjoint. What is important is that 

LjU = O, (5.3) 

where the right-hand side is the S x s zero matrix. Condition (j5.3p is, in fact, the solvability 
condition of the second equation in (|5.2p , and is equivalent to the condition that u/j satisfy 
the orthogonality relation ()2.1ip : 

(Z^,u^,) = 0. (5.4) 

Property (j5.4p is easily verified by substituting (j5.ip into the conservation law dQ/dt = 0. 

By ()5.4p . (j5.3p one can invert the second equation in (|5.2p and substitute the result in 
the first equation: 

LuR = -X^LJ^UR. (5.5) 

Recall that this generalized eigenvalue problem is considered on the restricted space (15. 4p . 
Both operators in (15. 5p are self-adjoint. Then, if L^^^ (and hence L/) is negative definite, 
then by Sylvester's law of inertia, the sign of is determined by the signs of the eigenvalues 
of L. If L is negative definite (again — on the restricted space (|5.4p . or equivalently, (j2.1ip ). 
then < and hence the solitary wave is dynamically linearly stable (see ()5.ip ). 

Thus, to summarize: The conditions of convergence of the ITEM and modified CGM 
coincide with the conditions under which the solitary wave is dynamically linearly stable 
if and only if operator L/ is negative definite. In particular, this occurs when u is a 
ground state. (Unfortunately, the latter fact is not readily determined by inspection for 
multi-component solitary waves.) Let us note that this result about the coincidence of 
convergence and stability conditions for ground-state solitary waves fully agrees with the 
results of [U [TO] , where it was proven that an ITEM- like method converges to ground states 
of (jl.2p and its generalization (j3.20p describing dynamics of Bose-Einstein condensates. 

Finally, we give a counterpart of the above statement for carrier wave equations using 
the KdV equation 

ut + 2uux + Uxxx = (5.6) 
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as an example. Its solitary wave u = u{x — ct) = u{^) with velocity c satisfies an equation 
— cu + Uxx = 0, whose linearized operator is L = 2u — c -\- d^- The iterative methods 
seeking a solitary wave with a prescribed value of power (jl.ip will converge provided that 
L is negative definite on a space of functions ^ satisfying a variant of (|2.1ip : 

(7z,V)=0. (5.7) 

On the other hand, the stability analysis of (j5.6p via an ansatz u{x, t) = u{^)+u{£,)e^'' leads 
to the eigenvalue problem dxLu = Xu. Upon the substitution w = d~^u, this eigenvalue 
problem is rewritten in the same form as ()5.2p |20) : 

Lu = Xw, Lw = —Xu, L = —dxLdx ■ (5-8) 

Note that due to the translational invariance of the solitary wave, u is a solution of Lu = 0. 
Thus, all considerations for the envelope equations carry over to the case of ()5.6p . and we 
conclude that the convergence conditions of the ITEM and CGM coincide with the stability 
conditions of the solitary wave if u is the ground state of L (or, equivalently, L is negative 
definite on the space defined by (|5.7p ). 

It may also be pointed out that in coupled multi-component generalizations of the KdV, 
there is only one parameter — the wave's velocity c — that is the analog of the propagation 
constant vector /2 for the envelope equations, like in (|2.ip or (|3.20p . Therefore, in this case, 
there is only one quadratic conserved quantity (which is probably the sum of the powers of 
all the components). In other words, s = 1, and hence the ITEM and modified CGM can 
converge only if p(L) < 1. 

6 Conclusions 

In this work, we obtained the convergence conditions of the iterative methods (12. 9p and (jA.ip 
that find multi-component solitary waves with prescribed values of quadratic conserved 
quantities (j2.6p . These convergence conditions are given by (j3.ip (provided that (jl.lOp holds 
for the ITEM (|2.9p ). which extend the convergence conditions of the single-component ITEM 
obtained in [T]. These conditions also turn out to be the same as the dynamical stability 
conditions for the nodeless (e.g., ground-state) solitary wave of the system of incoherently 
coupled NLS-type equations (j3.20p . which were obtained in [12J. For a single-component 
NLS-type equation (jl.2p . such a coincidence was observed in [1]. Earlier, similar statements 
were proven for ()1.2p and (I3.20p in [U |19] by different techniques. In Section 5 we showed 
that, in general, the convergence conditions of the iterative methods (j2.9p and (|A.ip . on 
one hand, and the dynamical stability conditions, on the other, coincide for ground-state 
solitary waves of all Hamiltonian nonlinear wave equations. 

Let us conclude with two remarks. First, even though we stated the ITEM (j2.9p in the 
main text of the paper while stating the CGM (jA.ip in Appendix, we remind the reader 
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(see Section 2) that if the ITEM converges slowly, then the modified CGM will provide 
considerable acceleration of the iterations. Alternatively, one can use the slowest mode 
elimination technique [2Tj to accelerate the algorithm of the ITEM. Comparison of these 
three methods was done in [11], with the modified CGM being found the fastest. 

Second, above we have explicitly mentioned the form of the operators employed by the 
ITEM and modified CGM for envelope equations (like (j2.ip and ()3.20p ) and for the carrier 
waves (like the KdV, (j5.6p ). For systems that couple envelope and carrier waves, which 
are commonly referred to as short-long wave interaction, or generalized Zakharov-Benney, 
equations, the formalism remains the same. 

Appendix: Modified CGM for solitary waves with a pre- 
scribed Q 

The steps of this algorithm for Eq. (j2.8ap are (operator C is defined in (j2.10p ): 

ro = N-^L^uo, do = ro - Uo {Uo, Uo)'^ {Uq, tq) , (A.la) 

(r„, Nd„) 



(k) (k) 
Un+1 = Un + andn , li^+i = U^+i 



\ Mk) ' K-l,...,SS^, 

(A.lc) 

Vn+i = N-1l(°)u„+i (A.ld) 

= K::7d;o ' ^^-'^^ 

dn+l = r„+i + Pndn — Un+1 (Z^n+li ^n+l)^"^ {^n+l, T^n+1 + Pndn) ■ (A. If) 

Equation (jA.laP defines the initial residual fq and the search direction dg. The first 
equation in (lA.lcP updates the iterative solution along the search direction d„ by making 
a "step" of "length" a„ found in (jATbj). Equations (jA.ll i.f) update the residual and the 
search direction using an auxiliary parameter /3„ computed in (jA.lep . 
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